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ABSTRACT 

This paper reports on the results of a numerical investigation designed to address how the initially anisotropic 
appearance of a GRB remnant is modified by the character of the circumburst medium and by the possible pres- 
ence of an accompanying supernova (SN). Axisymmetric hydrodynamical calculations of light, impulsive jets 
propagating in both uniform and inhomogeneous external media are presented, which show that the resulting 
dynamics of their remnants since the onset of the non-relativistic phase is different from the standard self- 
similar solutions. Because massive star progenitors are expected to have their close-in surroundings modified 
by the progenitor winds, we consider both free winds and shocked winds as possible external media for GRB 
remnant evolution. Abundant confirmation is provided here of the important notion that the morphology and 
visibility of GRB remnants are determined largely by their circumstellar environments. For this reason, their 
detectability is highly biased in favor of those with massive star progenitors; although, in this class of models, 
the beamed component may be difficult to identify because the GRB ejecta is eventually swept up by the ac- 
companying SN. The number density of asymmetric GRB remnants in the local Universe could be, however, 
far larger if they expand in a tenuous interstellar medium, as expected for some short GRB progenitor mod- 
els. In these sources, the late size of the observable, asymmetric remnant could extend over a wide, possibly 
resolvable angle and may be easier to constrain directly. 

Subject headings: gamma-rays: bursts - supernova remnants - hydrodynamics - shock waves - ISM: structure 



1. INTRODUCTION 

Relativistic jets are common in the astrophysical environ- 
ment. Objects known or sus pected to produce them include 
radio gala xies and quasars (Begelm an et al.|[l984| ), micro- 
quasars ([Mirabel & Ro driguez 1999 ) and gamma-ray bursts 
( Gorosabel et al . |2006j > . An important difference between jets 
of gamma-ray bursts (GRBs) and the better studied radio jets 
of quasars or microquasars is that active quasars often inject 
energy over extended periods of time into the jet while GRB 
sources are impulsive. Although quasar jets remain highly 
collimated throughout their lifetimes, GRB jets decelerate and 
expand significantly once they become nonrelativistic. Ex- 
pansion into a uniform medium has been well studied ( |Ayal| 
& Piran 2001), but the interaction of a GRB remnant with a 
nonuniform medium remains poorly understood. 

Much of our effort in this paper is therefore dedicated to 
determining how the morphology and dynamics of young 
GRB remnants is modified by the character of the circum- 
burst medium. Some of the questions at the forefront of atten- 
tion include the effects of the external medium and the degree 
to which GRB remnant dynamics and structures are modified 
by the presence of an accompanying supernova. We address 
both of these issues here. Because massive stars are expected 
to have their close-in surroundings modified by the progeni- 
tor winds, we consider both free winds and shocked winds as 
possible surrounding media for the GRB remnant evolution. 
Detailed hydrodynamic simulations of this interaction are pre- 
sented in §^4]and|6j while a brief description of the numerical 
methods and the initial models is giving in ^2] For complete- 
ness, the interaction with a constant-density medium is dis- 
cussed in ^3] The role of supernova explosions in shaping 
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the evolution and morphology of GRB remnants is discussed 
in ^5] The effects of a nonspherical circumburst medium are 
briefly addressed in Discussion and conclusions are pre- 
sented in ^7] 

2. NUMERICAL METHODS AND INITIAL MODEL 

2.1. The Underlying Dynamics 

For simplicity, lets consider a uniform GRB jet with sharp 
edges and a half-opening angle 6q, with an initial value of 6q. 
The typical angular size, R±, of the jet at t < to, where to is 
the time at which the jet's Lorentz factor 7 drops to 0$ 1 , is the 
same as for a spherical flow: 

R ± (t)ocl-^j J^ocf-J t^ k) tw (1) 

for an external den sity profile p ext = Ar~ k ( [Granot, Ramirez- 1 
Ruiz & Loeb |2 005 ). Here E is the true kinetic energy content 
of the jet, and £i S0 = f£ l E is the isotropic equivalent energy 
where ft ~ 9^/2 is the beaming factor. At these early times, 
the flow is described by the Blandford-McKee (1976) self- 
similar solution, which provides an accurate expression for 
the temporal evolution of the observed image size: 

/? ± (f) = 4xl0 16 (^y /8 (^y /8 cm (2) 

for constant density medium (k = 0) with n = 1 no cm -3 (p = 
1.67 x 10~ 24 po g cm -3 ), and 

/ 77 \ 1/4 / f x 3/4 

tf ± (0 = 2.5xl0 16 f^J rf^) cm (3) 

for a stellar wind environment (k = 2) with A*=A/5xl0 n g 
cm -1 . 

Notice that even in the most optimistic cases, the charac- 
teristic size of a GRB image is only of order 1 p as about 
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a day after the GRB at the Hubble distance (Gran ot & Loeb 
2001), and so it cannot be resolved by existing telescopes. 
Obviously, the challenge is made easier for nearby sources, 
where the late size of the observable remnant could extend 
over a wide, possibly resolvable angle (|Oren et al. |2004||Gra- 
not & Ram irez-Ruiz 2004 ; |Granot, Ramirez -Ruiz & Loeb 
2005). At late times t > %r, where %r is the nonrelativis- 



tic transition time, the jet is expected to gradually approach 
the Sedov-Taylor self-similar solution, asymptotically reach- 
ing R ± (t) oc (E/A) l ^ 5 - k) t 2 ^ 5 ~ k \ Att e <t< t NR there is, how- 
ever, a large uncertainty in the hydrodynamic evolution of the 
jet, and in particular in its rate of sideways expansion. 

To illustrate the importance of sideways expansion on the 
evolution of the observed image, lets consider two extreme 
assumption s: (i) relativistic lateral expan sion in the comov- 
ing frame ( |Rhoads|[T999[ |Sari et al.|^999] ) for which e « 
max(#o,7 -1 ) so that at to < t < t^ R we have 7 « 6q 1 « 
Oq 1 exp(-R/Ro), and (ii) little or no lateral expansion, Oq w Oq 
for t < %r, in which case appreciable lateral expansion oc- 
curs only when the jet becomes sub-relativistic and gradually 
approaches spherical symmetry. 

For relativistic lateral expansion, Oq ~ 1 at % R = t NR (E), 
where R NR (E) = ct^ R (E) = [(3 -k)E /AnAc 2 ] 1 '^ and the jet 
radius will be similar to that of the Sedov-Taylor solution, 
R ST (E,t) (Et 2 /A) l / (5 ~ k) . In this case, one expects the flow 
to approach spherical symmetry only after a few dynamical 
times. This is probably not the case as clearly illustrated by 
the results of numerical simulat ions ( [Granot et al.|2001||Can-| 
|nizzo et aL]|2004[ |Granot|2007| ) showing only modest lateral 
expansion as long as the jet is relativistic. If lateral expansion 
is neglected, the jet becomes sub-relativistic only at 



^NRffiso) - <^NR(£iso) = 



4ttAc 2 



l/(3-k) 



(4) 



For expansion in a constant density medium (k = 0), eq. ^ 
can be rewritten as 



^NR(^iso) - 0.3 



F x 1/3 / r-l X 1/3 
^51 \ I Jb 



P0 



30 



pc, 



while for a wind medium (k = 2), it becomes 



(Esi 
\A* 



30 



pc. 



(5) 



(6) 



From eq. ^ it follows that R^ R (E{ S0 ) is a factor of 

~ (E iso /E)^-« = f- l/(3 ~ k) ~ e- 2/Q - k) larger than R NR (E) = 

ct^ R (E) and a factor of ~ f^ l ^ 5 ~ k) ~ 0q 2 ^ 5 ~® larger than 
7?sT[£,%R(£iso)L Thus eq. simply states that the jet keeps 
its original opening angle, Vq « Oq until ^nr(^so), and hence 
at this time the jet is still far from being spherical. Thus, once 
the jet becomes sub-relativistic, we expect it to expand side- 
ways significantly, and become roughly spherical only when 
it has increased its radius by a factor of i/j. This should occur 
roughly at a time f sph when R ST (E,t sph ) = ipR NR (E iso ): 



%R 

-1/2 



(7) 



This is a factor of ~ w 14((9o/0.1) _1 larger than the ex- 
pected transition time for relativistic lateral expansion in the 
comoving frame. In the sections that follows we present quan- 
titative estimates of f sp h for GRB jets expanding in variety of 
circumburst environments. 
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Fig. 1. — Schematic plot illustrating the numerical initiation of the GRB 
explosion. 



2.2. Initial GRB Model 

Common to all calculations is the initiation of the GRB ex- 
plosion as two identical blobs expanding in opposite direc- 
tions into the circumburst medium. Calculations were done in 
two dimensions in cylindrical geometry using the PPM adap- 
tive mesh refinement code FLASH (ver 2.4). Both blobs and 
the circumburst medium are modeled by a cold, 7 = 5/3, ideal 
gas. The initial configuration is as follows. The computa- 
tional domain, as illustrated in Fig. [T] is a unprolonged cylin- 
der in which the ejecta move along the symmetry axis. In 
the inner region of each of the pancakes, the ejecta mass, Mj 
is distributed uniformly. In all runs vj = 0.3c, ARj/Rj ^0.5, 
<9j w 0.5, Mj w and R } w R NR (E lS0 ). 

Without a detailed understanding of the exact shape and en- 
ergy distribution of the ejecta, we have only an approximate 
description of how to constru ct the initial conditio ns. How- 
ever, as clearly illustrated by |Ayal & Piran| ( |2001| ), the late 
time evolution of the ejecta is rather insensitive to uncertain- 
ties in the initial conditions. We have considered various ini- 
tial densities, angular widths, and shapes of the collimated 
ejecta and found that these are indeed unimportant in deter- 
mining the late morphology of the remnant. This stems from 
the fact that at late times the mass of the remnant is dominated 
by the circumburst gas, which washes out any variations in the 
initial conditions of the ejecta. 

3. THE APPEARANCE OF A GRB REMNANT IN THE ISM 

In this section we present a quantitative discussion of how 
the GRB remnant morphology is modified by expansion into a 
uniform medium. Expansion into a constant density medium 
is expected in a variety of p rogenitor models, in particular 
those related to short GRBs ([Fryer et al.||1999{ |Bloom et al. 



[T9991 |Belczynski et al.||2006[ |Lee & Ramirez^Ru iz 2007). 
The discussion largely follows that of |Ayal & Piran| ( |2001| ), al- 
though analytical solutions are derived here to illustrate what 
may not be obvious from earlier derivations. 

In the absence of characteristic scales in stellar ejecta and in 
the ambient medium, self-similar, spherically symmetric solu- 
tions exist, and they are wid ely used to interpret ob servational 
data of GRB remn ants (e.g. Waxman et al . 1 1 998] ) . However, 
as argued in §2.1| near the non-relativistic transition the rem- 
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nant is far from being spherical and thus the Sedov-Taylor 
solution fails to provide an accurate description of the evo- 
lution of the GRB ejecta. A simple estimate for r sp h can be, 
however, obtain by assuming that the expansion of each of 
the collimated blobs is accurately described by a self-similar, 
spherically symmetric solution. Under these conditions, the 
problem reduces to depositing a finite amount of energy E at 
two different locations separated by a distance R } ~ R NR (E[ S0 ). 
This is valid as long as the external medium is uniform around 
the explosion site. The evolution of the shock radius for each 
of the blobs will then follow: 

/ E \ 1/5 

#ST(0 = d-J t 2 ' 5 (8) 

with ^=1.17 for 7 = 5/3. The ratio between the remnant 
width and height, 

, = *jgg> (9) 

«S T (0+*j' 

will approach unity as the two blobs expand and merge. The 
GRB remnant will then become nearly spherical in shape after 
a time 

(tf (€f fef- « 

when s = 0.9^0.9- Obviously, the above calculation is only 
sketchy and should be taken as an order of magnitude figure, 
as it not only assumes that the mass of the collimated ejecta is 
negligible with respect to the swept-up gas but also neglects 
the presence of shocked material throughout the interaction 
region. 

Detailed hydrodynamic simulations of the evolution of a 
GRB remnant in a uniform medium are presented in Fig. [2] 
where the density contours of the expanding collimated ejecta 
at various times in its hydrodynamical evolution are plotted. 
As the collimated ejecta collides with the ambient medium, 
a bow shock forms. The shock propagates in the direction 
of motion but also perpendicular to it and, over time, wraps 
around the expanding ejecta. Although initially the remnant 
may be highly nonspherical, the ratio between its width and 
height will approach unity as the two blobs expand, merge 
into a single structure and then finally become spherical in 
shape. The resultant structure will not be perfectly spherical 
as some density inhomogeneities around the equator result- 
ing from the encounter remain visible. It will be, however, 
difficult to distinguish it (based on morphology alone) from a 
supernova remnant after about 

(E \ 1/3 

4 ph ^3xl0 3 f -^J yr, (11) 

when ^ ~ 0.9. The governing parameters of the late evolu- 
tion of the remnant are the initial energ y of the jet, E, and 
the density of the ambient medium, p ( |Ayal & Piran 2001| ). 
As illustrated in Fig. |3j these two initial parameters can also 
determine the early evolution of the remnant for a fixed fy. 

4. EVOLUTION IN A CIRCUMSTELLAR WIND MEDIUM 

If the progenitors of GRBs are massive stars then there is 
an analogy to the explosions of core collapse supernovae, for 
which there is abundant evidence that they interact with the 
winds from the progenitor stars. In most supernova cases, 
the radial range that is observed is only out to a few pc, 
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Fig. 2. — The evolution of a GRB remnant in a constant density medium. 
The ejecta, £j = 10 50 erg, and surrounding ISM, po = 10 3 , are characterized by 
a 5/3 adiabatic index. Shown are logarithmic density cuts in g cm -3 . Calcu- 
lations were done in two-dimensional cylindrical coordinates for seven levels 
of refinement. The size of the computational domain was (0.3 pc) 2 . 




Fig. 3. — The evolution of a GRB remnant for p = 10 3 (left) and po = 1. 
Evolutionary ages in years are indicated in each frame together with corre- 
sponding size scales. As expected for a constant density medium, a unique 
combination of E, p and t has the dimensions of R. 
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Fig. 4. — The evolution of a GRB remnant in a 1/r 2 medium. The ejecta, 
£j = 10 51 erg, and surrounding stellar wind medium (v^ = 10 3 km s _1 and 

M = 2.5 x 1O -5 M0 yr -1 ) are characterized by a 5/3 adiabatic index. Shown 
is the evolution of the specific energy, e, in erg/g. Calculations were done 
in two-dimensional cylindrical coordinates for ten levels of refinement. The 
size of the computational domain was (16 pc) 2 . 

such that the mass loss characteristics have not changed sig- 
nificantly during the tim e that mass is supplied to the wind 
(Chevalie r & Li| |2000). The density in the wind depends 
on the type of progenitor. Red supergiant stars, which are 
thought to be the progenitors of Type II supernovae, have 
slow dense winds. Wolf-Rayet stars, which are believed to be 
the proge nitors of Type Ib/c supernovae and possibly o f long 
GRBs ( [Woosley|[T993l |MacFadyen & Woosley|^999] ), have 
faster, lower-density winds. The winds from WRs are char- 
acterized by mas s-loss rates M & 10~ 5 Mr., yr -1 and velocities 
v w « 10 3 km s -1 ( [Chiosi & Maeder|l986| . In a steady, spheri- 
cally symmetric wind, the stellar density is p QXt =Ar~ 2 , where 
A = 5 x lO n (M/lO- 5 M yr l )(v w /l0 3 km s" 1 )" 1 g cm" 3 . Note 
that for this choice of stellar wind parameters A* = 1. 

For this discussion we shall first assume the stellar wind 
is effectively spherical. The evolution of a GRB remnant in 
a stellar wind since the onset of the non-relativistic phase is 
summarized in Fig. [4] Similar resulting structures to those 
described in ^3] are clearly seen. A bow shock forms as each 
blob collides with the stellar wind, which eventually wraps 
around the ejecta before the two expanding shells collide to 
form a single structure. However, because in a wind medium 
the swept-up mass increases only linearly with radius, the 
GRB remnant decelerates much more slowly than in a uni- 
form medium. Moreover, in a wind medium, resistance to 
sideways expansion is increased. This is because the bow 
shock, as it wraps around the ejecta, encounters a steadily 
increasing ambient pressure. As a result, the remnant will 
become roughly spherical only after a time 



*sph : 



: 5 X 10 4 



E51 



(12) 



when $ ~ 0.9. Beyond this point, the evolution will evolve 
into a classical Sedov-Taylor supernova remnant evolution. 

The estimate given by equation ( [T2| ) could be inaccurate 
for a number of reasons. Depending upon the wind history 
of the progenitor star and the properties of the surrounding 




Fig. 5. — The evolution of a GRB remnant in a non-spherical 1/r 2 
medium. The ejecta, £j = 10 51 erg, and surrounding stellar wind medium 

(v esc = 10 3 km s _1 , C = 1 ,<f = 0.5, M = 2.5 x KT 5 M yr 1 , and Q = 0.8) are 
characterized by a 5/3 adiabatic index. Shown is the evolution of the specific 
energy, e, in erg/g. Calculations were done in two-dimensional cylindrical 
coordinates for seven levels of refinement. The size of the computational 
domain was (30 pc) 2 . 



ISM, the density structure around R sp h ~ 30pc could be quite 
complicated. The non-steady nature of the winds in massive 
stars together with the relatively large ISM pressure expected 
in star-forming regions, leaves open the possibility of interac- 
tion with denser material at much early times. In this case, 
the GRB remnant will start being decelerated by the exter- 
nal medium at a smaller radius than it would expanding into 
a free 1/r 2 wind. Much of our effort in q6| will therefore be 
dedicated to determining the contribution or the presupernova 
ejecta of Wolf-Rayet stars to the circumburst environment, 
and describing how this external matter can affect the observ- 
able characteristics of GRB remnants. 

Large-scale density gradients in the ambient medium could 
result in asymmetric, nonradial distortions. In this case, f sp h 
could be larger than that given by equation ( 12 ). For example, 
Fig. [5] shows the evolution of a GRB remnant in an asymmet- 
ric wind, under the assumption that the progenitor star expe- 
riences non-spherical mass loss close to critical rotation: in 
other words, a scenario in which a slower and denser wind 
is confined to the equatorial plane. To compute the latitudi- 
nal dependence of the wind properties of a star close to criti- 
cal rotation ideally requires multi-dimensional models of the 
star and its ou tflowing atmosphere, which are not available. 
|Langer| ([1998 ), however, argued that the stellar flux and the 
radius might still vary only weakly from pole to the equator 
in very luminous stars. We therefore applied equations simi- 
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lar to those found by Bjorkman & CassineTli| ([T993 ) for winds 
of rotating stars in the limit of large distance from the star: 



log 8 



3 (0) = Cv esc (l-ftsin0) v 



(13) 



where we set the parameters defined in Bjorkman & Cassinelli 
( |1993| ) to C = 1 W = 0-5, ft = Vrot/Vcrit = 0.8, and v crit = 
Vesc/V^ = [GM*(1 - k)R*)] 1 / 2 , with M* and R* being mass 
and radius of the star, and k standing for the ratio L/L Edd of 
stellar to Eddington luminosity. Under the above conditions, 
the GRB remnant expands more quickly and easily into the 
lower density wind at the poles, producing an increasingly 
asymmetric double-lobed structure. 

Finally, the estimate given by equation ( [T2| ) would be mod- 
ified if the beamed GRB is accompanied by an underlying 
supernova, as expected in the collapsar model (MacFadyen & 
Woosley|[T999l |MacFadyen et al. ||2Q0T| |Woosley & Bloom 
2006 ). The large-angle SN outflow, responsible for exploding 
the star and producing the 56 Ni, would generally carry more 
energy and inertia than the relativistic jet itself (e.g., [Sober 



bergetal. |2004| |Mazzali et aLpOOSl |Kaneko et al.|2007) , 
so that the latter always overtakes it, and sweeps up the GRB 
ejecta. It is to this problem that we now turn our attention. 

5. INTERACTION WITH AN UNDERLYING SUPERNOVAE 

It seems likely that GRBs originate in a very small fraction 
of massive stars that undergo a catastrophic energy release 
event toward the end of their evolution. Expressly, th e asso- 
ciation of some GRBs with type Ic supernovae (e.g.,[H j orth 
|et al. |[2003l |Stanek et al, ||2003l |Pian et al.||2006» has 
pointed a finger at deaths of massive Wolf -Ray et stars as the 
cause of GRBs, or at least a subset thereof. The central en- 
gine is b elieved to give rise to a polar outflow with two com- 
ponents (M acFadyen & Woosley 1999; Rami rez-Ruiz et aL] 
2002; Wo osley & Bloom|2006T rOne large-angle outflow (the 
SN), containing most of the energy and mass, is responsible 
for exploding the star and producing the 56 Ni to make the 
SN bright. A second outflow component (the GRB jet) occu- 
pies a narrower solid angle and probably contains less energy 
(which can range from comparable to much less), and most of 
its energy is in material with relativistic velocities (where the 
typical Lorentz factor of the material that carries most of the 
energy in this component can vary significantly among SN- 
GRBs). 

The large-angle SN outflow, carrying more energy and in- 
ertia than the relativistic jet itself, will generally sweep up the 
GRB ejecta before it has been much decelerated. An order 
of magnitude estimate for t SV h can be obtain by assuming that 
the dynamics of the laterally-expanding, GRB ejecta is ac- 
curately described by a self- similar, spherically solution in a 
1 /r 2 medium and that the SN outflow does not appreciably 
slow down. The evolution of the shock radius of the beamed 
remnant is then given by 

1/3 



#ST(0 = £w 



A/3 



(14) 



with £ w = 0.73 for 7 = 5/3. The GRB ejecta, as it clears 
the surrounding stellar matter, will be overtaking by the large 
scale-outflow SN at 

*" M0 (ir)(§?f *• 

where (3 = vsn /c<1. After this time, the merged system will 
quickly become spherical. Yet it is clear that this simple esti- 
mate is inadequate as a model for the real complex dynamics, 
which necessitates the use of hydrodynamical calculations. 
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Fig. 6. — The evolution of a GRB remnant interacting with an underlying 
(slower expanding) spherical supernova. The stellar wind parameters as well 
as the GRB ejecta initial quantities are the same as in Fig. [4] Shown is the 
evolution of the specific energy, e, in erg/g for two different SN explosion 
energies: 5 x 10 51 erg (left) and 5 x 10 50 erg (right). Calculations were done 
in two-dimensional cylindrical coordinates for ten levels of refinement. The 
size of the computational domain was (30 pc) 2 . 



The evolution of a GRB remnant accompanied by an un- 
derlying supernova is shown in Fig. [6] for two different SN 
explosion energies. The stellar wind parameters are chosen 
so that they are equal to those displayed in Fig. [4] Com- 
pared to Fig. [4] significant structural differences appear when 
the SN outflow has grown significantly in size and it starts to 
overtake the laterally expanding GRB ejecta. Two illustrative 
cases are depicted: £"sn- 5 x 10 51 erg (/^panels) and 5 x 10 50 
erg (right panels). In both cases, the SN outflow carries more 
momentum than the beamed ejecta and drives a blast wave 
that eventually sweeps up all the GRB-shocked medium. This 
happens before the two beamed blobs collide on the equa- 
tor. The merged system will become roughly spherical very 
soon after. For E$n = 5 x 10 51 erg (5 x 10 50 erg), we obtain 
hph ~ 250 yr (3 x 10 3 yr). Not surprisingly, the presence of an 
underlying spherical supernova seriously modifies the simple 
estimate given by equation ( [T2| ) and limits our ability to de- 
cipher the presence of a beamed component in a GRB explo- 
sion. 

Prevailing to all these calculations is the initiation of the 
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Fig. 7. — The pressure structure, along the z axis of a GRB remnant inter- 
acting with an underlying spherical supernova. The simulation is the same as 
in Fig.|6]for % = 5 x 10 50 erg and t = 51, 190, 580 yr. 



explosion as a pressure-driven blast wave by deposition of the 
explosion energy, £sn, entirely as thermal energy. In the inner 
region, an ejecta mass, Msn ~ 5M , are distributed uniformly. 
This may accurately model the Sedov-Taylor stage of SN rem- 
nant evolution after the ratio of swept-up mass t o the mass of 
the or iginal stellar ejecta exceeds roughly 19 ( [Fabian et al. 
|1983| ). In most cases, deceleration of the SN outflow will be- 
gin only sometime after the GRB ejecta has been swept up. As 
clearly seen in Fig. [7] the SN outflow at this early stage is not 
accurately described by the Sedov-Taylor solution. However, 
numerous tests show that our results are not strongly depen- 
dent upon the assumed mass of the SN ejecta or on whether 
kinetic energy rather than thermal energy is distributed (such 
that the velocity profile is linear; similar to the Sedov solution) 
in the inner region. None of these complications is likely to 
seriously modified our estimate for t sp h- 



6. GRBs INSIDE PRE-EXISITING WIND-DRIVEN BUBBLES 

So far we have considered either the uniform ambient 
medium case or the 1 /r 2 wind case on its own. However, 
since the winds in massive stars are non- steady, the density 
structure is more complex ( Wijers 2001, Chevali er et al. 
|2004l|Ramrrez-Ruiz et al.|2005t|van Marie et al.|2006| ). The 
preburst stellar wind depends on the ev olutionary stages prior 
to (an d during) the Wolf-Rayet stage ( |Ramirez-Ruiz et a~ 
|2001| ). For Galactic stars, a standard evolutionary track is 
to start as an O star, evolve through a red supergiant (RSG) 
phase or luminous blue variable (LBV) phase with consid- 
erable mass loss, and ending as a Wolf-Rayet star ( [Garcia- 
jSegura et aL]|1996a[ |Garcia-Segura et al. || 1996b)). At low 



metallicity, the RSG phase may be absent; this may also be 
the case for some binary stars. 

As an example, we follow the dynamics of a GRB remnant 
around a 35 M star (as calculated by |Chevalier et al. |2004| ), 
which evolves (at solar metallicity) through a long-lived RSG 
stage with prominent consequences for the evolution of the 
circumstellar matter. The wind velocity in the Wolf-Rayet 
phase is 10 3 km s -1 , and the mass-loss rate 1O 5 M yr -1 . The 
ISM pressure and density are assumed to be typical of the 
hot, low-density phase of a starburst galaxy, with P{ sm ~ 10 7 
K cm -3 and a density of 4 x 10~ 25 g cm -3 . When the fast 
Wolf-Rayet wind v w ~ 10 3 km s -1 starts blowing, it sweeps 
up the RSG wind material into a shell. The termination shock 
of the Wolf-Rayet wind is located at R t w 0.4 pc and RSG 
shell at R YSg ^1.7 pc. Because the pressure in the shocked 
wind is nearly in equilibrium with the ISM, and the temper- 




FiG. 8. — The evolution of a GRB remnant inside the wind bubble structure 
expected around a 35 Mq massive star. The GRB ejecta initial quantities are 
the same as in Fig. [4] Shown is the evolution of the specific energy, e, in erg/g 
at t = 51,660,2301) and 3400 yr. The individual frames have been succes- 
sively rotated by tt/2. Calculations were done in two-dimensional cylindri- 
cal coordinates for eight levels of refinement. The size of the computational 
domain was (10 pc) 2 . 




radius (cm) 



Fig. 9. — The density structure, along the z axis of a GRB remnant inter- 
acting with a wind bubble structure. The simulation is the same as in Fig. [8] 
for t = 5 1 , 660, 2300 and 3400 yr. 



ature ~ lO^v^/lO 3 km s l ) 2 K, the density in the bubble is 
- 8 x 10 25 (Pi sm /10 7 )(v w /10 3 kms" 1 )- 2 g cm" 3 , independent of 
the mass-loss rate and the ambient density. The extent of the 
constant density region is ~ 4R t out to the dense red super- 
giant shell. 

The resulting evolution a beamed GRB remnant in the cir- 
cumstellar medium expected around a 35 M massive stel- 
lar progenitor is summarized in Fig. [8] The presence of the 
sharp density gradient will only affect the dynamics of the 
GRB remnant when it size is comparable to, or exceeds, the 
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scale length of the gradient. Before this time, the evolution 
of the remnant is similar to that depicted in Fig. [4] The 
bow shock propagates in the direction of motion but also per- 
pendicular to it and, over time, wraps around the expanding 
ejecta. In a wind bubble, the shock front will expand within 
the stellar wind until it reaches the sharp density discontinu- 
ity at about 1.7 pc. The encounter with the RSG shell happens 
before the GRB remnant has time to expand laterally, which 
allows for an elongation of the RSG shell in the z axis. A 
less pronounced elongation is also seen perpendicular to the z 
axis, which results from the collision on the equator of the two 
beamed blobs. Superimposed on this large-scale deformation 
one can also notice the familiar effects caused by the devel- 
opment of Rayleigh-Taylor (R-T) instability. The instability 
grows rather quickly and within the following 10 3 yr several 
elongated spikes extend from the shell. At w 2 x 10 3 yr, the 
remnant has a mean radius of 3 pc (Fig. [9]) and the spikes 
are even more pronounced than before. Both shell and R-T 
spikes advance with a velocity of 300-500 km s" 1 , which is 
a very small fraction of the random velocities observed in the 
hot cavity of the remnant. These fast motions are induced by 
inflection of consecutive shock waves ramming the irregular 
shell. The pressure in the cavity is almost uniform, whereas 
the density varies chaotically. The knots of ejecta formed by 
the R-T instability are believed to be responsible for several 
forms of observable radiation in some young SNRs, includ- 
ing radio synchrotron and optical emission lines. The sheared 
motions resulting from the instability could lead to an ampli- 
fication of the magnetic field strengt h and enhance the bright 
radio remnant ( Jun & Norm anp996| ). 

The growth time of the instability may be roughly estimated 
in the following way. The interaction of the GRB remnant 
with the wind cavity leads to an increase of the shell-driving 
pressure by 



AP = 



2nRL 



The mass of the shell is equal to 

4 



^rsg = -7T# 3 sg p rS g. 



(16) 



(17) 



AP causes an acceleration g rsg of the shell, satisfying the 
equation 



47rtf 2 A/> = M rsg g rsg . 



(18) 



With no other forces involved, the contact surface separating 
both fluids is unstable to perturbations of all wavelengths and 
the fluids interpenetrate. The instability grows exponentially 
on a characteristic time scale 



'4^Atf^\ 1/2 
3£i 



(19) 



where A is the perturbation wavelength, corresponding to a 
growth time of 



fo-T ; 



10 3 



2 pc 



A 



0.2 pc 



1/2 



o l/2 

rsg,- 



r 1/2 

21^51 y r > 



(20) 



where E51 is the total energy of the GRB remnant and 
Prsg =10~ 2 Vrsg,-2i gem -3 . Indeed, the observed spikes have 
dimensions of tenths of parsecs and grow on a time scale com- 
parable to the estimated one. For this discussion we have as- 
sumed that the remnant evolution is effectively adiabatic and 
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Fig. 10. — The evolution of a GRB remnant accompanied by an underlying 
spherical supernova. The wind bubble structure as well as the GRB ejecta 
initial quantities are the same as in Fig. [8] Shown is the evolution of the 
specific energy, e, in erg/g at t = 380,5707765 and 950 yr. The individual 
frames have been successively rotated by tt/2. Calculations were done in 
two-dimensional cylindrical coordinates for eight levels of refinement. The 
size of the computational domain was (20 pc) 2 . 



should be modified to include the effects of radiative cool- 
ing 4 , which are expected to become important after about few 
thousand years. 

The disruption of the RSG shell is, as expected, very sen- 
sitive to the total amount of energy and momentum released 
by the GRB explosion. To illustrate this, we study the evo- 
lution of a GRB remnant accompanied by an underlying 10 52 
erg SN in a wind bubble structure such as that illustrated in 
Fig. [8] The large-angle SN outflow, carrying significantly 
more energy and inertia than the relativistic jet itself, sweeps 
up the GRB ejecta before it is decelerated by the RSG shell. 
The resultant shell, clearly apparent in the first frame of Fig. 
[TO] taken at £=380 yr, will then be pushed outward at faster 
velocity than it would be in the absence of the SN ejecta. 
In subsequent evolutionary phases one can observe a broad- 
ening of the merged shell, accompanied by a decreased of 
its density, and the development of small- amplitude density 
perturbations. After t = 3180 yr, the remnant has becomes 
highly irregular already when it has grown to w 13 pc, and 
is also rather weakly elongated in the direction of motion of 
the beamed ejecta (Fig. [TT]). There are also multiple kine- 
matic components within the remnant. Fast-moving knots and 
fast-moving flocculi, dense fragments of SN ejecta, expand- 
ing from the explosion center with velocities of several thou- 
sand km s" 1 . Much slower (several hundred km s -1 ) flocculi 
are clearly shocked and accelerated remnants of stellar mate- 
rial ejected by the SN progenitor prior to the explosion. 

4 For example, the growth of the R-T ins tability in the radiating s hell is 
found to be higher than in the adiabatic case ( Chevalier & Blondi njl995} . 
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Fig. 11. — The density (p), specific energy (e), velocity (v), and Mach 
number (M) of a GRB remnant accompanied by an underlying spherical su- 
pernova (Msn ~ interacting with a wind bubble structure. The simu- 
lation is the same as in Fig.QoJbut for t = 3180 yr. 

The calculations above demonstrate how the measurable 
properties of GRB remnants depend on the nature of the pro- 
genitor star and the medium around it. Moreover, counts of 
GRB remnants as a function of age may have huge selec- 
tion effects, as the actual age of the GRB may be consider- 
ably less than the kinematical age estimated from the radius 
of the filaments divided by the expansion velocity. Although 
nurture makes a huge difference, hydrodynamical models of 
young GRB remnants are also sensitive to the structure of the 
ejecta, which might be rather complicated in detail as demon- 
strate d by inspection of numeri cal models of SN explosions 
(e.g. |Chevalier & Blondin[T995] ). 



7. GRB REMNANTS AND THEIR DETECTABILITY 
7.1. Asymmetric GRB Remnants in the Local Universe 

It is obvious from the discussions in this paper that the dy- 
namics of GRB remnants are complex, especially because 
of rich interactions between the ejecta and the circumburst 
medium. The structure of ejecta is also important, in par- 
ticular for young GRB remnants. Abundant confirmation was 
provided of the important notion that the morphology and vis- 
ibility of GRB remnants are determined largely by their cir- 
cumstellar environment, from the initial density gradient cre- 
ated by the progenitor to the effects of large- and small-scale 
circumburst structures in later evolution. This distribution is 
expected to be nonuniform around massive GRB progenitors 
because of the significant mass loss and the dynamical effects 
of the stellar winds. 

The presence of a density gradient will only affect the dy- 
namics of a GRB when the remnant size is comparable to, or 



TABLE 1 

Mean space density of asymmetric GRB remnants. 



^sph 



0s P h[3£grb = IGpc 3 yr l ] 



k = 
k = 2 



fvrsg ' 



' 1 pc 



E„ 



^,51 



3xl034/V /3 yr 
5 x lO^A; 1 yr 
lO 2 ^; 1 yr 
3 x 10 3 yr 
< 10 3 yr 



3xl0 3 4 /3 po 1/3 A" 1 %bGpc- 3 
5 x 10 4 E 5l A~ l f-% rh Gpc- 3 



10 2 £ 5 iA;7^grb Gpc- 3 
3 x lO 3 /;^ Gpc" 3 
<10 3 /- 1 $R grb Gpc- 3 



exceeds, the scale length of the gradient. Before this time, 
the density can be treated as approximately uniform. In the 
absence of characteristic scales in stellar ejecta and in the am- 
bient medium, self- similar, spherically symmetric solutions 
exist, and they are widely used to interpret observational data 
on young GRB remnants. However, even for the simplest den- 
sity distributions, we found that the resulting structure and dy- 
namics of their remnants are very different from the standard 
self- similar solutions. This is mainly because at early stages 
the morphology of a beamed GRB remnant would be very 
different from that of a spherical explosion. In principle, this 
can be used to identify those remnants although the dynami- 
cal complexity of their sourrounding circumburst medium se- 
riously limits ourr ability to decipher their presence, in partic- 
ular around massive star progenitors. 

The values of the mean space density, sp h, of asymmetric 
GRB remnants expanding into a variety of circumburst envi 
ronments are given in Table 1 . Taking the long GRB rate ex 
pected in the local universe to be 3? gr b = 0.3 Gpc -3 yr -1 ([Guetta 
& Piran 2007), an upper limit to the number density 01 asym 
metric remnants arising from massive stellar progenitors is 



^sph^xlO-^siA^^Mpc- 



(21) 



This is a generous upper limit, since it assumes that all long 
GRBs occur in a free 1/r 2 stellar wind and that the relativistic 
component is energetically dominant. For comparison, an ul- 
traconservative lower limit to the space density of long GRBs 
gives 



^sph - 



10_5 io5 Mpc " 3 - 



(22) 



With the rates given by Guetta & Piran ( 2005), the frequency 
of short GRBs is about 20 times higher than that of long 
GRBs. The much lower external density (likely uniform) 
medium expected around short GRBs (e.g., |Lee et al.|[2 005) 
suggests 



■ 6x 10" J £ 5 f p. 



3r7 l/3 -1/3 fb 1 



10 



Mpc 



(23) 



where p-3 = 10 3 po. Thus the rate of asymmetric short GRB 



remnants could be much higher than equation ( 22 ). However, 
the remnant's visibility would be highly biased in favor of 
those with massive progenitors. 



7.2. Observational Prospects 

The most difficult task at present is to relate hydrodynami- 
cal modeling to observations. A few of the observables, such 
as expansion rates and thicknesses of the flow structures, can 
be relatively easily determined from the models. However, 
modeling radio and X-ray emission is in general difficult, as 
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Fig. 12. — Compilation of GRB, supernova and SN remnant X-ray light 
curves (0.3- 10 keV) pres ented as (isotropic) luminosity distances as a func - 
tion of age (adapted from Kouveliotou et al. 2004 Immler & Kuntz 2005). 
Superimposed on these plots, on a common scale, is the schematic lightcurves 
of two GRB remnants expanding into different external density environments. 
The inset figures show the evolution of the GRB remnant in a 1/r 2 density 
profile (the simulation is the same as in Fig. [4j. The grey curve shows the 
evolution of a GRB remnant in a pre-existing wind-bubble (the simulation is 
the same as in Fig. [8j. To estimate the X-ray lightcurves, we have summed 
up the internal energy density of the gas, e- mt , for each zone in the simulations 
and plotted Lx ~ fx^int as a function of time, where fx ~ 0.01 is the fraction 
of internal energy that is radiated away in the 0.3-10 keV range. 



we are still lacking an understanding of how electrons are ac- 
celerated in shocks. Very similar difficulties are encountered 
in modeling nonthermal X-rays. Thermal X-ray spectra are 
in principle easier to model, but in practice the difficulties 
are formidable. The reason for these difficulties is our poor 
understanding of a number of topics, such as the amount of 
electron heating in collisionless shocks, the detailed structure 
and composition of ejecta, their clumping, the presence of 
the inhomogeneous circumstellar medium, and the presence 
of dust. 

The difficult task of interpreting observations with the help 
of hydrodynamical models is perhaps best illustrated by Cas 
A. This is a remnant of a massive star explosion and a 
classic prototype shell supernova remnant. It has bee n de- 
tected throughout the whole electromagnetic spectrum (|Ryle 
& SmithH 19481 |Ashworth|[T980l |Aharonian et aL||2001i |Fe- 
sen|2001||Hwang et al.|2004| ). Observations in various wave- 
length bands probe very different components of the remnant: 
synchrotron radio emission gives us information about rela- 
tivistic electrons, thermal X-ray emission is produced by the 
bulk of the shocked hot gas, much cooler gas in radiative 
shocks emits at optical wavelengths, and observations in in- 
frared reveal still cooler gas and dust. However, we still do not 
understand what is the relationships bet ween all these featu res 
and the remnant's hydrodynamics (e.g. Hwang et al. 2004| ). 

To estimate the emissivity of GRB remnant without under- 
taking the complicated effort of calculating X-ray spectra, we 
have summed up the internal energy density of the gas, 6i n t, 



for each zone in the simulations and plotted L x 
a function of time in Fig. 12 where fx ~ 0.01 



~ f x e int as 
is the frac- 
tion of internal energy that is radiated away in the 0.3-10 keV 
range. This is only a very approximate procedure and should 
be taken as an order of magnitude estimate at present. Super- 
imposed on these plots, on a common scale, are all GRB after- 
glows with X-ray luminosity measurements covering several 
tens of days. There are, unfortunately, only a few curves, be- 
cause such measurements can only be made on GRBs that are 
relatively nearby, but their light curves should be illustrative. 
We then compared these X-ray lightcurves with those of su- 
pernovae and historical Galactic SN remnants. The resulting 
plot is striking in several ways. Despite the huge disparity in 
initial appearance, there are indications of a common conver- 
gence of all classes of phenomena to a common resting place: 
L x ~ 10 39 - 10 40 erg s -1 about a few years after the explosive 
event. Moreover, it clearly illustrates that the transition from 
a GRB to a SN remnant appears to be rather smooth. Clearly, 
detailed studies relating hydrodynamical modeling to obser- 
vations are needed to study the transition from a GRB into a 
stellar remnant. 

Fig. [T2| shows the schematic lightcurves of two GRB rem- 
nants expanding into different external density environments. 
For al/r 2 density profile, as in the black curve of Fig. 12 the 



luminosity of the remnant is initially dominated by the emis- 
sion of the individual beamed components as they interact 
with the stellar wind. The luminosity continues with a quasi- 
steady decay rate until the individual beamed components col- 
lide to form a single structure. The resultant lightcurve will 
then be characterized by a modest increase in luminosity. For 
a GRB expanding inside a pre-existing wind-bubble, as in the 
grey curve of Fig. [12| the resultant lightcurves ca n evolve 
much faster into luminous remnants such as Cas A (Hwang 
et al.|2004| ) or W49B ( jMiceli et aLp0 06 ) due to their strong 
interaction with the dense circumburst medium. 

This plot summarized many of the issues outlined in this 
paper, in which we argued that the morphology and visibility 
of GRB remnants are governed by the pre-existing structure 
of their birthplace environments. Since GRB remnants result 
from the impact of their ejecta with circumstellar gas, their 
visibility is highly biased in favor of those with massive pro- 
genitors. Many young GRBs from massive progenitors would 
be bright because their ejected mass is interacting with nearby 
gas expelled by the progenitor itself. This circumstellar gas is 
likely to have mass comparable to that of the accompanied SN 
debris and will not extend much further than a few parsecs. 
After less than a century, the blast wave from the GRB will 
pass through this relatively dense circumstellar gas. Inferring 
the the presence of a beamed component in these hypernova 
scenarios would be challenging. The number density of asym- 
metric GRB remnants in the local Universe could be far larger 
if they expand in a tenuous interstellar medium as expected, 
for example, in the merger of two neutron stars (although 
there are reasons to suspect that the ejecta may not be too 
narrowly beamed) and may be easier to constrain directly (ac- 
knowledging the obvious trade-offs in sensitivity and angular 
resolution, particularly for radio and X-ray observations). 



We have benefited from many useful discussions with C. 
Fryer, J. Granot, T. Piran and S. Woosley. We are especially 
grateful to W. Zhang for countless insightful conversations. 
The software used in this work was in part developed by 
the DOE-supported ASCI/Alliance Center for Astrophysical 
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